knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE 

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix = paste0('_final_with_plan_and_barriers', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.1649 0.0000 5.0000 275

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan', # todo: do we need this??
    'barriers_self_cw', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'barriers_self_cb', 
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily individual's experienced persuasion",  
    "Daily partner's experienced persuasion", 
    "Daily individual's experienced pressure", 
    "Daily partner's experienced pressure", 
    "Daily individual's experienced pushing", 
    "Daily partner's experienced pushing", 
    "Day", 
    "Own actionplan",
    'Partner actionplan',
    "Daily barriers",
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean individual's experienced persuasion", 
    "Mean partner's experienced persuasion", 
    "Mean individual's experienced pressure", 
    "Mean partner's experienced pressure", 
    "Mean individual's experienced pushing", 
    "Mean partner's experienced pushing", 
    'Mean barriers',
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily individual's experienced persuasion)", 
  "sd(Daily partner's experienced persuasion)", # OR partner received
  "sd(Daily individual's experienced pressure)", 
  "sd(Daily partner's experienced pressure)", 
  "sd(Daily individual's experienced pushing)", 
  "sd(Daily partner's experienced pushing)", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,12),
  "Between-Person Effects" = c(13,20),
  "Random Effects" = c(21, 27), 
  "Additional Parameters" = c(28,28)
  )


rows_to_pack_ordinal <- list(
  "Within-Person Effects" = c(2+5,12+5),
  "Between-Person Effects" = c(13+5,20+5),
  "Random Effects" = c(21+5, 27+5), 
  "Additional Parameters" = c(28+5,28+6)
  )

Self-Reported MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | dd | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | dd | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  control = list(adapt_delta = 0.90),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
if (check_models) {
  check_brms(pa_sub, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_sub, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.34 [1.29, 1.39]         1.16      0.75
##  persuasion_partner_cw 1.08 [1.06, 1.13]         1.04      0.92
##       pressure_self_cw 1.05 [1.02, 1.09]         1.02      0.96
##    pressure_partner_cw 1.06 [1.04, 1.10]         1.03      0.94
##        pushing_self_cw 1.03 [1.01, 1.08]         1.02      0.97
##     pushing_partner_cw 1.14 [1.11, 1.18]         1.07      0.88
##     persuasion_self_cb 1.04 [1.02, 1.09]         1.02      0.96
##  persuasion_partner_cb 3.65 [3.48, 3.84]         1.91      0.27
##       pressure_self_cb 3.62 [3.45, 3.81]         1.90      0.28
##    pressure_partner_cb 2.34 [2.24, 2.45]         1.53      0.43
##        pushing_self_cb 2.33 [2.23, 2.44]         1.53      0.43
##     pushing_partner_cb 3.29 [3.13, 3.45]         1.81      0.30
##       barriers_self_cw 3.31 [3.15, 3.47]         1.82      0.30
##       barriers_self_cb 1.01 [1.00, 1.55]         1.00      0.99
##                    day 1.13 [1.10, 1.18]         1.06      0.88
##  Tolerance 95% CI
##      [0.72, 0.77]
##      [0.89, 0.95]
##      [0.92, 0.98]
##      [0.91, 0.97]
##      [0.92, 0.99]
##      [0.84, 0.90]
##      [0.92, 0.98]
##      [0.26, 0.29]
##      [0.26, 0.29]
##      [0.41, 0.45]
##      [0.41, 0.45]
##      [0.29, 0.32]
##      [0.29, 0.32]
##      [0.64, 1.00]
##      [0.85, 0.91]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         91%
##     lognormal          9%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##  neg. binomial (zero-infl.)         84%
##               beta-binomial          9%
##                   lognormal          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 3 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 5, observations = 1776, p-value = 0.28
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.003378378
## sample estimates:
## outlier frequency (expected: 0.00154842342342342 ) 
##                                        0.002815315
if (do_priorsense) {
  priorsense_vars <- c(
      'Intercept',
      'b_persuasion_self_cw',
      'b_persuasion_partner_cw',
      'b_pressure_self_cw',
      'b_pressure_partner_cw',
      'b_pushing_self_cw',
      'b_pushing_partner_cw'
  )
  
  hurdle_priorsense_vars <- c(
    'Intercept_hu',
    'b_hu_persuasion_self_cw',
    'b_hu_persuasion_partner_cw',
    'b_hu_pressure_self_cw',
    'b_hu_pressure_partner_cw',
    'b_hu_pushing_self_cw',
    'b_hu_pushing_partner_cw'
  )
  
  gc()
  priorsense::powerscale_sensitivity(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_dens(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_ecdf(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
  priorsense::powerscale_plot_quantities(pa_sub, variable = c(priorsense_vars, hurdle_priorsense_vars))
}
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report, rope_range
## = rope_range_continuous, : Coefficients were exponentiated. Double check if
## this was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
Intercept 2.49*** 0.60 [ 1.53, 4.04] 1.000 [0.84, 1.20] 0.002 32.035 Very Strong Evidence 1.000 33697 28937 49.28*** 3.96 [41.95, 57.78] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 12605 22320
Within-Person Effects
Daily individual’s experienced persuasion 1.47*** 0.14 [ 1.23, 1.79] 1.000 [0.84, 1.20] 0.012 >100 Overwhelming Evidence 1.000 53988 30427 1.04 0.03 [ 0.99, 1.10] 0.926 [0.93, 1.08] 0.892 0.032 Strong Evidence for Null 1.000 30055 31124
Daily partner’s experienced persuasion 1.32** 0.13 [ 1.10, 1.63] 0.999 [0.84, 1.20] 0.148 4.490 Moderate Evidence 1.000 45996 29424 1.02 0.02 [ 0.98, 1.07] 0.831 [0.93, 1.08] 0.983 0.014 Very Strong Evidence for Null 1.000 42578 33714
Daily individual’s experienced pressure 0.80 0.20 [ 0.47, 1.47] 0.806 [0.84, 1.20] 0.361 0.202 Moderate Evidence for Null 1.000 37908 25268 0.89 0.05 [ 0.78, 1.00] 0.971 [0.93, 1.08] 0.238 0.068 Strong Evidence for Null 1.000 43179 31403
Daily partner’s experienced pressure 1.49 0.57 [ 0.81, 4.82] 0.892 [0.84, 1.20] 0.232 0.306 Weak Evidence for Null 1.000 24485 20576 0.95 0.05 [ 0.86, 1.05] 0.846 [0.93, 1.08] 0.708 0.014 Very Strong Evidence for Null 1.000 45159 33810
Daily individual’s experienced pushing 1.16 0.17 [ 0.87, 1.57] 0.851 [0.84, 1.20] 0.571 0.123 Moderate Evidence for Null 1.000 41350 30381 0.99 0.03 [ 0.93, 1.05] 0.669 [0.93, 1.08] 0.972 0.009 Very Strong Evidence for Null 1.000 41853 37229
Daily partner’s experienced pushing 1.58** 0.26 [ 1.17, 2.45] 0.998 [0.84, 1.20] 0.034 6.451 Moderate Evidence 1.000 32558 24454 0.98 0.03 [ 0.92, 1.04] 0.737 [0.93, 1.08] 0.965 0.010 Very Strong Evidence for Null 1.000 45057 33647
Day 0.97 0.21 [ 0.63, 1.49] 0.561 [0.84, 1.20] 0.588 0.110 Moderate Evidence for Null 1.000 77129 29862 1.01 0.07 [ 0.89, 1.16] 0.566 [0.93, 1.08] 0.725 0.008 Very Strong Evidence for Null 1.000 77457 30134
Own actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.21*** 0.06 [ 1.10, 1.33] 1.000 [0.84, 1.20] 0.398 34.893 Very Strong Evidence 1.000 87289 28688 1.04** 0.01 [ 1.01, 1.07] 0.998 [0.93, 1.08] 0.993 0.518 Weak Evidence for Null 1.000 80040 28881
Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.87 0.47 [ 0.30, 2.51] 0.602 [0.84, 1.20] 0.255 0.281 Moderate Evidence for Null 1.000 23084 29360 1.06 0.17 [ 0.77, 1.47] 0.652 [0.93, 1.08] 0.340 0.016 Very Strong Evidence for Null 1.000 22190 28337
Mean partner’s experienced persuasion 1.26 0.67 [ 0.43, 3.60] 0.667 [0.84, 1.20] 0.238 0.297 Moderate Evidence for Null 1.000 23496 28628 0.98 0.16 [ 0.71, 1.34] 0.561 [0.93, 1.08] 0.364 0.016 Very Strong Evidence for Null 1.000 22153 28525
Mean individual’s experienced pressure 0.30 0.22 [ 0.06, 1.19] 0.956 [0.84, 1.20] 0.049 1.564 Weak Evidence 1.000 33009 31559 0.66 0.26 [ 0.30, 1.47] 0.852 [0.93, 1.08] 0.087 0.042 Strong Evidence for Null 1.000 9042 15981
Mean partner’s experienced pressure 0.25 0.18 [ 0.05, 1.01] 0.974 [0.84, 1.20] 0.032 2.267 Weak Evidence 1.000 35470 30743 0.52 0.21 [ 0.23, 1.15] 0.950 [0.93, 1.08] 0.039 0.093 Strong Evidence for Null 1.000 8878 15559
Mean individual’s experienced pushing 1.13 0.90 [ 0.23, 5.53] 0.563 [0.84, 1.20] 0.176 0.408 Weak Evidence for Null 1.000 30160 30199 1.55 0.45 [ 0.86, 2.77] 0.931 [0.93, 1.08] 0.071 0.052 Strong Evidence for Null 1.000 12963 21788
Mean partner’s experienced pushing 1.49 1.22 [ 0.30, 7.30] 0.686 [0.84, 1.20] 0.155 0.460 Weak Evidence for Null 1.000 30808 30057 1.57 0.46 [ 0.87, 2.81] 0.938 [0.93, 1.08] 0.062 0.058 Strong Evidence for Null 1.000 12584 21865
Mean barriers 1.35* 0.16 [ 1.07, 1.72] 0.994 [0.84, 1.20] 0.152 1.511 Weak Evidence 1.000 72775 30710 1.08* 0.04 [ 1.01, 1.16] 0.991 [0.93, 1.08] 0.460 0.190 Moderate Evidence for Null 1.000 46784 33848
Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.20 0.18 [0.91, 1.63] NA NA NA NA NA 1.000 24322 29367 0.30 0.04 [0.22, 0.40] NA NA NA NA NA 1.000 17652 23878
sd(Daily individual’s experienced persuasion) 0.16 0.13 [0.01, 0.44] NA NA NA NA NA 1.000 17188 21318 0.11 0.03 [0.06, 0.17] NA NA NA NA NA 1.000 21505 25370
sd(Daily partner’s experienced persuasion) 0.22 0.14 [0.02, 0.51] NA NA NA NA NA 1.000 15551 18313 0.06 0.03 [0.01, 0.12] NA NA NA NA NA 1.000 16566 15276
sd(Daily individual’s experienced pressure) 0.43 0.40 [0.02, 1.79] NA NA NA NA NA 1.000 14632 21652 0.07 0.06 [0.00, 0.26] NA NA NA NA NA 1.000 22600 22413
sd(Daily partner’s experienced pressure) 0.77 0.59 [0.04, 2.34] NA NA NA NA NA 1.000 14057 19351 0.06 0.05 [0.00, 0.19] NA NA NA NA NA 1.000 22879 23239
sd(Daily individual’s experienced pushing) 0.42 0.23 [0.04, 0.94] NA NA NA NA NA 1.000 12577 14194 0.05 0.04 [0.00, 0.14] NA NA NA NA NA 1.000 16328 20593
sd(Daily partner’s experienced pushing) 0.42 0.29 [0.03, 1.15] NA NA NA NA NA 1.000 11434 16306 0.05 0.04 [0.00, 0.14] NA NA NA NA NA 1.000 17218 21743
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 0.67 0.01 [0.65, 0.70] NA NA NA NA NA 1.000 64881 30448 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.75), b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.71), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.8). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Hurdle part of the model on the left, non-zero part towards the right side of the table

conds_eff <- conditional_spaghetti(
  pa_sub, 
  effects = c(
    'persuasion_self_cw',
    'persuasion_partner_cw',
    'pressure_self_cw',
    'pressure_partner_cw',
    'pushing_self_cw',
    'pushing_partner_cw'
  ),
  x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  ),
  group_var = 'coupleID',
  plot_full_range = TRUE,
  y_limits = c(0, 100),
  y_label = "Same-Day MVPA",
  y_labels = c('Probability of Being Active', 'Minutes of MVPA When Active', 'Overall Expected Minutes of MVPA'),
  , filter_quantiles = .9995
  , font_family = 'Candara'
)
## This is posterior version 1.6.0
## 
## Attaching package: 'posterior'
## The following object is masked from 'package:bayesplot':
## 
##     rhat
## The following objects are masked from 'package:stats':
## 
##     mad, sd, var
## The following objects are masked from 'package:base':
## 
##     %in%, match
## Registering fonts with R
## Scanning ttf files in C:\Windows/Fonts, C:\Users\pascku\AppData\Local/Microsoft/Windows/Fonts ...
## Extracting .afm files from .ttf files...
## C:\Windows\Fonts\Candara.ttf : Candara already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarab.ttf : Candara-Bold already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarai.ttf : Candara-Italic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaral.ttf : Candara-Light already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candarali.ttf : Candara-LightItalic already registered in fonts database. Skipping.
## C:\Windows\Fonts\Candaraz.ttf : Candara-BoldItalic already registered in fonts database. Skipping.
## Found FontName for 0 fonts.
## Scanning afm files in C:/Users/pascku/AppData/Local/R/cache/R/renv/cache/v5/windows/R-4.4/x86_64-w64-mingw32/extrafontdb/1.0/a861555ddec7451c653b40e713166c6f/extrafontdb/metrics
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
## Warning: Dropping 'draws_df' class as required metadata was removed.
print(conds_eff)

$persuasion_self_cw

## Warning: Removed 144 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 228 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00776

$persuasion_partner_cw

## Warning: Removed 94 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 199 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.00688

$pressure_self_cw

## Warning: Removed 52 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 91 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0124

$pressure_partner_cw

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Picking joint bandwidth of 0.0277

$pushing_self_cw

## Picking joint bandwidth of 0.009

$pushing_partner_cw

## Picking joint bandwidth of 0.0128

Note. This graphic illustrates the relationship between social control and moderate to vigorous physical activity (MVPA) using a Bayesian Hurdle-Lognormal Multilevel Model. The predictor is centered within individuals to examine how deviations from their average social control relate to same-day MVPA. Shaded areas indicate credible intervals, thick lines show fixed effects, and thin lines represent random effects, highlighting variability across couples. The plots display the probability of being active, expected minutes of MVPA when active, and combined predicted MVPA. The bottom density plot visualizes the posterior distributions of slope estimates, transformed to represent multiplicative changes in odds ratios (hurdle component) or expected values. Medians and 95% credible intervals (2.5th and 97.5th percentiles) are shown. Effects are significant, when the 95% credible interval does not overlap 1.

x_label = c(
    'Received Persuasion',
    'Exerted Persuasion',
    'Received Pressure',
    'Exerted Pressure',
    'Received Plan-Related Pushing',
    'Exerted Plan-Related Pushing'
  )

home_dir <- getwd()
output_dir <- file.path(home_dir, 'Output', 'Plots')

for (i in 1:length(conds_eff)) {
  effname <- names(conds_eff)[i]
  eff_plot <- conds_eff[[i]]
  x_label_i <- x_label[[i]]
  
  rmarkdown::render(
    file.path(output_dir, 'BeautifulPlotWithNote.Rmd'), 
    output_file = file.path(output_dir, paste0('Graphic_', effname, '.pdf')),
    params = list(
      home_dir = home_dir,
      output_dir = output_dir,
      p_i = eff_plot,
      p_name = effname,
      x_label = x_label_i
      ),
    envir = new.env(),
    quiet = TRUE
  )
}

print('done')

Comparing effect size of pressure and pushing

hypothesis(pa_sub, "pressure_self_cw < pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... < 0    -0.11      0.08    -0.23     0.02      11.53
##   Post.Prob Star
## 1      0.92     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner + 
    barriers_self_cw + barriers_self_cb + 
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
if (check_models) {
  check_brms(pa_obj_log, log_pp_check = TRUE)
  DHARMa.check_brms.all(pa_obj_log, integer = TRUE, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.11 [1.07, 1.17]         1.05      0.90
##  persuasion_partner_cw 1.15 [1.11, 1.21]         1.07      0.87
##       pressure_self_cw 1.08 [1.05, 1.14]         1.04      0.93
##    pressure_partner_cw 1.10 [1.07, 1.16]         1.05      0.91
##        pushing_self_cw 1.17 [1.13, 1.23]         1.08      0.85
##     pushing_partner_cw 1.16 [1.12, 1.22]         1.08      0.86
##       pressure_self_cb 3.31 [3.10, 3.55]         1.82      0.30
##    pressure_partner_cb 3.51 [3.28, 3.76]         1.87      0.29
##       barriers_self_cw 1.01 [1.00, 1.37]         1.01      0.99
##       barriers_self_cb 1.07 [1.04, 1.14]         1.04      0.93
##                    day 1.06 [1.03, 1.13]         1.03      0.94
##       weartime_self_cw 1.06 [1.03, 1.12]         1.03      0.95
##       weartime_self_cb 1.21 [1.16, 1.27]         1.10      0.83
##  Tolerance 95% CI
##      [0.85, 0.93]
##      [0.83, 0.90]
##      [0.88, 0.96]
##      [0.86, 0.94]
##      [0.81, 0.89]
##      [0.82, 0.90]
##      [0.28, 0.32]
##      [0.27, 0.30]
##      [0.73, 1.00]
##      [0.88, 0.96]
##      [0.89, 0.97]
##      [0.89, 0.97]
##      [0.79, 0.86]
## 
## Moderate Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 5.97 [5.56, 6.42]         2.44      0.17
##  persuasion_partner_cb 6.44 [5.99, 6.92]         2.54      0.16
##        pushing_self_cb 5.19 [4.84, 5.58]         2.28      0.19
##     pushing_partner_cb 5.13 [4.78, 5.51]         2.26      0.20
##  Tolerance 95% CI
##      [0.16, 0.18]
##      [0.14, 0.17]
##      [0.18, 0.21]
##      [0.18, 0.21]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         88%
##     lognormal          9%
##       weibull          3%
## 
## Predicted Distribution of Response
## 
##                Distribution Probability
##                   lognormal         72%
##                     tweedie         16%
##  neg. binomial (zero-infl.)          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 2 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 14, observations = 1594, p-value < 2.2e-16
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0006273526 0.0062735257
## sample estimates:
## outlier frequency (expected: 0.00321831869510665 ) 
##                                        0.008782936
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(pa_obj_log, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(pa_obj_log, variable = priorsense_vars)
}
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report, :
## Coefficients were exponentiated. Double check if this was intended.
summary_pa_obj %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 123.18*** 7.40 [109.03, 139.10] 1.000 [0.94, 1.06] 0.000 >100 Overwhelming Evidence 1.000 10287 18834
Within-Person Effects
Daily individual’s experienced persuasion 1.02 0.02 [ 0.99, 1.06] 0.881 [0.94, 1.06] 0.987 0.014 Very Strong Evidence for Null 1.000 32647 32365
Daily partner’s experienced persuasion 1.03 0.02 [ 0.99, 1.07] 0.940 [0.94, 1.06] 0.970 0.023 Very Strong Evidence for Null 1.000 43952 31850
Daily individual’s experienced pressure 0.95 0.04 [ 0.86, 1.04] 0.870 [0.94, 1.06] 0.556 0.014 Very Strong Evidence for Null 1.000 43518 30052
Daily partner’s experienced pressure 0.97 0.04 [ 0.90, 1.04] 0.814 [0.94, 1.06] 0.763 0.010 Very Strong Evidence for Null 1.000 56254 32810
Daily individual’s experienced pushing 1.03 0.03 [ 0.97, 1.08] 0.818 [0.94, 1.06] 0.904 0.012 Very Strong Evidence for Null 1.000 42843 33705
Daily partner’s experienced pushing 1.01 0.02 [ 0.96, 1.06] 0.629 [0.94, 1.06] 0.982 0.007 Very Strong Evidence for Null 1.000 50585 31498
Day 0.97 0.05 [ 0.88, 1.07] 0.746 [0.94, 1.06] 0.693 0.007 Very Strong Evidence for Null 1.000 76155 29860
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.03** 0.01 [ 1.01, 1.05] 0.999 [0.94, 1.06] 0.999 0.457 Weak Evidence for Null 1.000 77752 27470
Daily weartime 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.06] 1.000 1.711 Weak Evidence 1.000 82611 30094
Between-Person Effects
Mean individual’s experienced persuasion 1.02 0.16 [ 0.74, 1.39] 0.547 [0.94, 1.06] 0.310 0.014 Very Strong Evidence for Null 1.000 9251 17057
Mean partner’s experienced persuasion 0.87 0.14 [ 0.63, 1.19] 0.814 [0.94, 1.06] 0.210 0.022 Very Strong Evidence for Null 1.000 9170 16781
Mean individual’s experienced pressure 1.05 0.19 [ 0.74, 1.49] 0.603 [0.94, 1.06] 0.261 0.010 Very Strong Evidence for Null 1.000 14364 23206
Mean partner’s experienced pressure 1.04 0.18 [ 0.74, 1.47] 0.595 [0.94, 1.06] 0.271 0.010 Very Strong Evidence for Null 1.000 12665 21665
Mean individual’s experienced pushing 1.12 0.26 [ 0.71, 1.78] 0.696 [0.94, 1.06] 0.190 0.014 Very Strong Evidence for Null 1.000 12416 20531
Mean partner’s experienced pushing 1.29 0.29 [ 0.83, 2.05] 0.874 [0.94, 1.06] 0.112 0.025 Very Strong Evidence for Null 1.000 12435 20072
Mean barriers 1.05 0.03 [ 0.99, 1.10] 0.958 [0.94, 1.06] 0.744 0.042 Strong Evidence for Null 1.000 45006 32464
Mean weartime 1.00* 0.00 [ 1.00, 1.00] 0.983 [0.94, 1.06] 1.000 0.090 Strong Evidence for Null 1.000 37030 33323
Random Effects
sd(Intercept) 0.32 0.04 [0.25, 0.43] NA NA NA NA NA 1.000 12078 20334
sd(Daily individual’s experienced persuasion) 0.05 0.02 [0.02, 0.09] NA NA NA NA NA 1.000 19548 14612
sd(Daily partner’s experienced persuasion) 0.05 0.02 [0.00, 0.10] NA NA NA NA NA 1.000 11945 12399
sd(Daily individual’s experienced pressure) 0.06 0.06 [0.00, 0.23] NA NA NA NA NA 1.000 16266 23012
sd(Daily partner’s experienced pressure) 0.04 0.03 [0.00, 0.13] NA NA NA NA NA 1.000 25114 22322
sd(Daily individual’s experienced pushing) 0.09 0.03 [0.02, 0.16] NA NA NA NA NA 1.000 12694 10874
sd(Daily partner’s experienced pushing) 0.05 0.03 [0.00, 0.12] NA NA NA NA NA 1.000 14159 17910
Additional Parameters
sigma 0.54 0.01 [0.52, 0.56] NA NA NA NA NA 1.000 57239 28922 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.82), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.75). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

range(df_double$aff, na.rm = T)
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
if (check_models) {
  check_brms(mood_gauss, log_pp_check = FALSE)
  DHARMa.check_brms.all(mood_gauss, integer = FALSE)
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.17 [1.13,  1.23]         1.08      0.85
##  persuasion_partner_cw 1.14 [1.10,  1.20]         1.07      0.88
##       pressure_self_cw 1.13 [1.09,  1.19]         1.06      0.88
##    pressure_partner_cw 1.07 [1.04,  1.13]         1.04      0.93
##        pushing_self_cw 1.26 [1.21,  1.33]         1.12      0.79
##     pushing_partner_cw 1.13 [1.09,  1.19]         1.07      0.88
##       pressure_self_cb 4.88 [4.55,  5.23]         2.21      0.21
##    pressure_partner_cb 4.77 [4.46,  5.11]         2.18      0.21
##       barriers_self_cw 1.01 [1.00,  1.84]         1.00      0.99
##       barriers_self_cb 1.07 [1.04,  1.13]         1.04      0.93
##                    day 1.07 [1.04,  1.13]         1.03      0.93
##  Tolerance 95% CI
##      [0.81, 0.89]
##      [0.84, 0.91]
##      [0.84, 0.92]
##      [0.88, 0.96]
##      [0.75, 0.83]
##      [0.84, 0.91]
##      [0.19, 0.22]
##      [0.20, 0.22]
##      [0.54, 1.00]
##      [0.88, 0.96]
##      [0.88, 0.96]
## 
## Moderate Correlation
## 
##                   Term  VIF    VIF 95% CI Increased SE Tolerance
##     persuasion_self_cb 9.80 [9.12, 10.55]         3.13      0.10
##  persuasion_partner_cb 9.64 [8.97, 10.37]         3.11      0.10
##        pushing_self_cb 7.06 [6.58,  7.58]         2.66      0.14
##     pushing_partner_cb 6.98 [6.51,  7.50]         2.64      0.14
##  Tolerance 95% CI
##      [0.09, 0.11]
##      [0.10, 0.11]
##      [0.13, 0.15]
##      [0.13, 0.15]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        cauchy         44%
##        normal         41%
##   exponential          6%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##        tweedie         41%
##  beta-binomial         38%
##    half-cauchy          6%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 14, observations = 1776, p-value =
## 2.096e-05
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.004316152 0.013190801
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.007882883
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(mood_gauss, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(mood_gauss, variable = priorsense_vars)
}
summary_mood <- summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Sampling priors, please wait...
summary_mood %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.80*** 0.11 [ 3.58, 4.02] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.001 8225 15991
Within-Person Effects
Daily individual’s experienced persuasion 0.01 0.02 [-0.04, 0.06] 0.620 [-0.11, 0.11] 1.000 0.005 Very Strong Evidence for Null 1.000 51887 32038
Daily partner’s experienced persuasion 0.01 0.02 [-0.04, 0.06] 0.678 [-0.11, 0.11] 1.000 0.005 Very Strong Evidence for Null 1.000 59637 31864
Daily individual’s experienced pressure -0.05 0.06 [-0.18, 0.08] 0.763 [-0.11, 0.11] 0.833 0.006 Very Strong Evidence for Null 1.000 45767 28568
Daily partner’s experienced pressure -0.06 0.06 [-0.19, 0.06] 0.844 [-0.11, 0.11] 0.791 0.009 Very Strong Evidence for Null 1.000 44947 28359
Daily individual’s experienced pushing 0.04 0.04 [-0.04, 0.11] 0.837 [-0.11, 0.11] 0.975 0.008 Very Strong Evidence for Null 1.000 49871 34302
Daily partner’s experienced pushing 0.11* 0.04 [ 0.02, 0.18] 0.993 [-0.11, 0.11] 0.542 0.156 Moderate Evidence for Null 1.000 41442 29456
Day 0.24** 0.08 [ 0.09, 0.39] 0.999 [-0.11, 0.11] 0.044 0.615 Weak Evidence for Null 1.000 72175 29192
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 0.13*** 0.02 [ 0.10, 0.16] 1.000 [-0.11, 0.11] 0.137 >100 Overwhelming Evidence 1.000 78022 28614
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.14 0.29 [-0.44, 0.72] 0.689 [-0.11, 0.11] 0.266 0.015 Very Strong Evidence for Null 1.000 7958 14644
Mean partner’s experienced persuasion 0.43 0.29 [-0.15, 1.00] 0.928 [-0.11, 0.11] 0.103 0.042 Strong Evidence for Null 1.000 7666 14449
Mean individual’s experienced pressure -0.18 0.30 [-0.78, 0.41] 0.728 [-0.11, 0.11] 0.242 0.011 Very Strong Evidence for Null 1.000 10073 19492
Mean partner’s experienced pressure -0.31 0.30 [-0.90, 0.28] 0.852 [-0.11, 0.11] 0.172 0.015 Very Strong Evidence for Null 1.000 9729 18338
Mean individual’s experienced pushing 0.28 0.41 [-0.54, 1.12] 0.753 [-0.11, 0.11] 0.170 0.015 Very Strong Evidence for Null 1.000 10182 17943
Mean partner’s experienced pushing -0.07 0.41 [-0.89, 0.77] 0.564 [-0.11, 0.11] 0.208 0.012 Very Strong Evidence for Null 1.000 10015 17318
Mean barriers 0.36*** 0.04 [ 0.28, 0.44] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 49900 32196
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.61 0.08 [0.48, 0.80] NA NA NA NA NA 1.000 11221 20168
sd(Daily individual’s experienced persuasion) 0.04 0.03 [0.00, 0.12] NA NA NA NA NA 1.000 13475 17798
sd(Daily partner’s experienced persuasion) 0.04 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 14426 19079
sd(Daily individual’s experienced pressure) 0.07 0.06 [0.00, 0.27] NA NA NA NA NA 1.000 21788 21202
sd(Daily partner’s experienced pressure) 0.09 0.07 [0.00, 0.28] NA NA NA NA NA 1.000 18175 20881
sd(Daily individual’s experienced pushing) 0.07 0.05 [0.00, 0.18] NA NA NA NA NA 1.000 13031 16597
sd(Daily partner’s experienced pushing) 0.10 0.05 [0.01, 0.21] NA NA NA NA NA 1.000 13954 13987
Additional Parameters
sigma 0.88 0.02 [0.85, 0.91] NA NA NA NA NA 1.000 62268 29718 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.9), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.76), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.73), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.73), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.75), b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.83). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 7) 

hist(log(df_double$reactance+0.1), breaks = 10)

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
if (check_models) {
  check_brms(reactance_ordinal)
  DHARMa.check_brms.all(reactance_ordinal, outliers_type = 'bootstrap')
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##       pressure_self_cw 3.88 [3.68, 4.10]         1.97      0.26
##    pressure_partner_cw 1.42 [1.37, 1.48]         1.19      0.70
##        pushing_self_cw 1.32 [1.27, 1.38]         1.15      0.76
##     pushing_partner_cw 1.16 [1.12, 1.21]         1.08      0.86
##     persuasion_self_cb 1.11 [1.07, 1.15]         1.05      0.90
##  persuasion_partner_cb 1.05 [1.03, 1.10]         1.03      0.95
##       pressure_self_cb 1.33 [1.29, 1.39]         1.15      0.75
##    pressure_partner_cb 1.17 [1.13, 1.22]         1.08      0.85
##        pushing_self_cb 4.49 [4.25, 4.75]         2.12      0.22
##       barriers_self_cw 4.50 [4.26, 4.76]         2.12      0.22
##       barriers_self_cb 4.80 [4.54, 5.08]         2.19      0.21
##                    day 2.08 [1.98, 2.18]         1.44      0.48
##  Tolerance 95% CI
##      [0.24, 0.27]
##      [0.67, 0.73]
##      [0.73, 0.78]
##      [0.83, 0.89]
##      [0.87, 0.93]
##      [0.91, 0.97]
##      [0.72, 0.78]
##      [0.82, 0.88]
##      [0.21, 0.24]
##      [0.21, 0.23]
##      [0.20, 0.22]
##      [0.46, 0.50]
## 
## Moderate Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 7.84 [7.40, 8.32]         2.80      0.13
##  persuasion_partner_cw 9.24 [8.71, 9.81]         3.04      0.11
##     pushing_partner_cb 5.73 [5.41, 6.07]         2.39      0.17
##  Tolerance 95% CI
##      [0.12, 0.14]
##      [0.10, 0.11]
##      [0.16, 0.18]
## Error in h(simpleError(msg, call)) : 
##   error in evaluating the argument 'x' in selecting a method for function 'print': Predictive errors are not defined for ordinal or categorical models.
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 7 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
##  DHARMa bootstrapped outlier test
## 
## data:  model.check
## outliers at both margin(s) = 0, observations = 486, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004115226
## sample estimates:
## outlier frequency (expected: 0.000473251028806584 ) 
##                                                   0
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(reactance_ordinal, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(reactance_ordinal, variable = priorsense_vars)
}
summary_reactance_ordinal <- summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) 
## Sampling priors, please wait...
summary_reactance_ordinal %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 3.82*** 1.39 [ 1.93, 8.02] 1.000 [0.84, 1.20] 0.001 48.727 Very Strong Evidence 1.000 23250 26841
Intercept[2] 8.22*** 3.08 [ 4.05, 17.81] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 23129 26259
Intercept[3] 22.70*** 9.21 [ 10.67, 52.42] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 23388 25422
Intercept[4] 95.42*** 45.43 [ 39.79, 257.89] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 24833 25331
Intercept[5] 1637.00*** 1393.12 [368.24, 11665.00] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 39723 30374
Within-Person Effects
Daily individual’s experienced persuasion 0.73** 0.09 [ 0.57, 0.92] 0.995 [0.84, 1.20] 0.120 2.389 Weak Evidence 1.000 25197 26372
Daily partner’s experienced persuasion 1.01 0.13 [ 0.77, 1.28] 0.534 [0.84, 1.20] 0.836 0.057 Strong Evidence for Null 1.000 34318 26552
Daily individual’s experienced pressure 1.76* 0.39 [ 1.07, 2.76] 0.985 [0.84, 1.20] 0.052 0.996 Weak Evidence for Null 1.000 22085 22183
Daily partner’s experienced pressure 1.23 0.40 [ 0.52, 2.44] 0.737 [0.84, 1.20] 0.326 0.098 Strong Evidence for Null 1.000 18106 13983
Daily individual’s experienced pushing 1.31* 0.16 [ 1.02, 1.69] 0.982 [0.84, 1.20] 0.240 0.553 Weak Evidence for Null 1.000 33937 29293
Daily partner’s experienced pushing 0.90 0.13 [ 0.66, 1.22] 0.754 [0.84, 1.20] 0.668 0.073 Strong Evidence for Null 1.000 31840 24011
Day 1.37 0.64 [ 0.53, 3.48] 0.746 [0.84, 1.20] 0.242 0.064 Strong Evidence for Null 1.000 40322 29394
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.00 0.10 [ 0.83, 1.21] 0.504 [0.84, 1.20] 0.935 0.051 Strong Evidence for Null 1.000 39235 30250
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.33 0.93 [ 0.33, 5.61] 0.660 [0.84, 1.20] 0.187 0.083 Strong Evidence for Null 1.000 16836 22501
Mean partner’s experienced persuasion 2.29 1.83 [ 0.51, 12.15] 0.854 [0.84, 1.20] 0.109 0.125 Moderate Evidence for Null 1.000 16531 21647
Mean individual’s experienced pressure 2.65 2.16 [ 0.53, 13.99] 0.884 [0.84, 1.20] 0.087 0.128 Moderate Evidence for Null 1.000 20515 28252
Mean partner’s experienced pressure 0.85 0.70 [ 0.15, 4.01] 0.582 [0.84, 1.20] 0.173 0.061 Strong Evidence for Null 1.000 17717 24538
Mean individual’s experienced pushing 0.76 0.77 [ 0.11, 6.13] 0.608 [0.84, 1.20] 0.137 0.073 Strong Evidence for Null 1.000 17231 22556
Mean partner’s experienced pushing 0.04* 0.06 [ 0.00, 0.48] 0.994 [0.84, 1.20] 0.005 1.934 Weak Evidence 1.000 18520 24475
Mean barriers 0.93 0.25 [ 0.55, 1.64] 0.610 [0.84, 1.20] 0.469 0.088 Strong Evidence for Null 1.000 20462 24232
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.08 0.28 [0.61, 1.74] NA NA NA NA NA 1.000 11347 19261
sd(Daily individual’s experienced persuasion) 0.33 0.19 [0.02, 0.71] NA NA NA NA NA 1.001 5828 9273
sd(Daily partner’s experienced persuasion) 0.18 0.14 [0.01, 0.53] NA NA NA NA NA 1.000 12913 14882
sd(Daily individual’s experienced pressure) 0.49 0.32 [0.03, 1.30] NA NA NA NA NA 1.000 9173 11304
sd(Daily partner’s experienced pressure) 0.66 0.53 [0.03, 2.05] NA NA NA NA NA 1.001 8324 12613
sd(Daily individual’s experienced pushing) 0.23 0.18 [0.01, 0.65] NA NA NA NA NA 1.000 8043 14635
sd(Daily partner’s experienced pushing) 0.17 0.15 [0.01, 0.63] NA NA NA NA NA 1.000 16092 18227
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.78), b_Intercept[4] and b_Intercept[3] (r = 0.84),
##   b_pressure_partner_cb and b_persuasion_partner_cb (r = 0.74). This might
##   lead to inappropriate results. See 'Details' in '?rope'.

conditional_spaghetti(
  reactance_ordinal, 
  effects = c('persuasion_self_cw', 'pressure_self_cw')
  , group_var = 'coupleID'
  #, n_groups = 15
  , plot_full_range = T
)

\(persuasion_self_cw <img src="01_FinalModelsBarriers-ONLY_WITH_PLAN-_files/figure-html/report_reactance_ordinal-3.png" width="2400" />\)pressure_self_cw

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    barriers_self_cw + barriers_self_cb + 
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
if (check_models) {
  check_brms(is_reactance)
  DHARMa.check_brms.all(is_reactance, integer = FALSE)
}
## # Check for Multicollinearity
## 
## Low Correlation
## 
##                   Term  VIF   VIF 95% CI Increased SE Tolerance
##     persuasion_self_cw 1.09 [1.06, 1.14]         1.04      0.92
##  persuasion_partner_cw 1.23 [1.18, 1.28]         1.11      0.82
##       pressure_self_cw 1.04 [1.02, 1.10]         1.02      0.96
##    pressure_partner_cw 1.04 [1.02, 1.10]         1.02      0.96
##        pushing_self_cw 1.11 [1.07, 1.15]         1.05      0.90
##     pushing_partner_cw 1.23 [1.19, 1.29]         1.11      0.81
##     persuasion_self_cb 2.96 [2.81, 3.13]         1.72      0.34
##  persuasion_partner_cb 3.42 [3.24, 3.61]         1.85      0.29
##       pressure_self_cb 2.09 [1.99, 2.20]         1.45      0.48
##    pressure_partner_cb 2.06 [1.97, 2.17]         1.44      0.48
##        pushing_self_cb 2.45 [2.33, 2.58]         1.56      0.41
##     pushing_partner_cb 2.29 [2.18, 2.41]         1.51      0.44
##       barriers_self_cw 1.06 [1.04, 1.11]         1.03      0.94
##       barriers_self_cb 1.17 [1.14, 1.22]         1.08      0.85
##                    day 1.05 [1.02, 1.10]         1.02      0.95
##  Tolerance 95% CI
##      [0.88, 0.94]
##      [0.78, 0.84]
##      [0.91, 0.98]
##      [0.91, 0.98]
##      [0.87, 0.93]
##      [0.78, 0.84]
##      [0.32, 0.36]
##      [0.28, 0.31]
##      [0.46, 0.50]
##      [0.46, 0.51]
##      [0.39, 0.43]
##      [0.41, 0.46]
##      [0.90, 0.97]
##      [0.82, 0.88]
##      [0.91, 0.98]
## # Distribution of Model Family
## 
## Predicted Distribution of Residuals
## 
##  Distribution Probability
##        normal         34%
##          beta         16%
##        cauchy         16%
## 
## Predicted Distribution of Response
## 
##   Distribution Probability
##      bernoulli         75%
##  beta-binomial         16%
##       binomial          9%
## 
## Divergences:
## 0 of 40000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 40000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 29 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.

## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  model.check
## outliers at both margin(s) = 1, observations = 486, p-value = 0.6217
## alternative hypothesis: true probability of success is not equal to 0.001998002
## 95 percent confidence interval:
##  0.0000520929 0.0114105115
## sample estimates:
## frequency of outliers (expected: 0.001998001998002 ) 
##                                          0.002057613
if (do_priorsense) {
  gc()
  priorsense::powerscale_sensitivity(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_dens(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_ecdf(is_reactance, variable = priorsense_vars)
  priorsense::powerscale_plot_quantities(is_reactance, variable = priorsense_vars)
}
summary_is_reactance <- summarize_brms(
  is_reactance, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.57 0.22 [0.26, 1.24] 0.923 [0.83, 1.20] 0.138 0.186 Moderate Evidence for Null 1.000 28387 31479
Within-Person Effects
Daily individual’s experienced persuasion 0.69* 0.11 [0.48, 0.92] 0.993 [0.83, 1.20] 0.093 1.997 Weak Evidence 1.000 23677 25761
Daily partner’s experienced persuasion 1.12 0.20 [0.79, 1.66] 0.740 [0.83, 1.20] 0.606 0.095 Strong Evidence for Null 1.000 27636 26371
Daily individual’s experienced pressure 2.31* 1.12 [1.03, 9.08] 0.978 [0.83, 1.20] 0.047 0.970 Weak Evidence for Null 1.000 15749 15577
Daily partner’s experienced pressure 1.79 1.26 [0.44, 13.34] 0.817 [0.83, 1.20] 0.146 0.230 Moderate Evidence for Null 1.000 18210 17447
Daily individual’s experienced pushing 1.43* 0.23 [1.06, 2.07] 0.990 [0.83, 1.20] 0.129 1.092 Weak Evidence 1.000 24749 27471
Daily partner’s experienced pushing 0.88 0.21 [0.54, 1.50] 0.695 [0.83, 1.20] 0.483 0.109 Moderate Evidence for Null 1.000 30479 27591
Day 1.29 0.69 [0.45, 3.69] 0.679 [0.83, 1.20] 0.236 0.065 Strong Evidence for Null 1.000 42497 31844
Own actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.03 0.12 [0.83, 1.29] 0.612 [0.83, 1.20] 0.878 0.063 Strong Evidence for Null 1.000 39766 29361
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 3.61 3.66 [0.51, 30.37] 0.904 [0.83, 1.20] 0.064 0.251 Moderate Evidence for Null 1.000 18329 25297
Mean partner’s experienced persuasion 6.05 6.76 [0.74, 67.12] 0.953 [0.83, 1.20] 0.034 0.409 Weak Evidence for Null 1.000 17821 25403
Mean individual’s experienced pressure 42.14* 78.46 [1.43, 2281.97] 0.985 [0.83, 1.20] 0.008 1.319 Weak Evidence 1.000 16155 22563
Mean partner’s experienced pressure 0.65 1.39 [0.01, 35.77] 0.580 [0.83, 1.20] 0.066 0.161 Moderate Evidence for Null 1.000 12462 22258
Mean individual’s experienced pushing 0.48 0.79 [0.02, 12.72] 0.674 [0.83, 1.20] 0.079 0.122 Moderate Evidence for Null 1.000 15358 21859
Mean partner’s experienced pushing 0.03 0.06 [0.00, 1.03] 0.974 [0.83, 1.20] 0.012 0.767 Weak Evidence for Null 1.000 15681 24453
Mean barriers 1.10 0.42 [0.54, 2.40] 0.598 [0.83, 1.20] 0.361 0.116 Moderate Evidence for Null 1.000 19221 25224
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.95 0.42 [1.25, 2.94] NA NA NA NA NA 1.000 12148 22042
sd(Daily individual’s experienced persuasion) 0.47 0.21 [0.06, 0.93] NA NA NA NA NA 1.001 7019 7320
sd(Daily partner’s experienced persuasion) 0.40 0.24 [0.03, 0.99] NA NA NA NA NA 1.000 11120 11475
sd(Daily individual’s experienced pressure) 1.20 0.68 [0.12, 2.93] NA NA NA NA NA 1.001 8521 8928
sd(Daily partner’s experienced pressure) 1.64 1.03 [0.15, 4.12] NA NA NA NA NA 1.001 13076 12068
sd(Daily individual’s experienced pushing) 0.31 0.23 [0.02, 0.86] NA NA NA NA NA 1.000 9740 15750
sd(Daily partner’s experienced pushing) 0.42 0.31 [0.02, 1.22] NA NA NA NA NA 1.000 13238 14532
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pressure_self_cw', 'pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

\(pressure_self_cw <img src="01_FinalModelsBarriers-ONLY_WITH_PLAN-_files/figure-html/report_is_reactance-3.png" width="2400" />\)pushing_self_cw 

hypothesis(is_reactance, "exp(pressure_self_cw) > exp(pushing_self_cw)")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (exp(pressure_sel... > 0     1.52      2.67    -0.44     5.37       5.31
##   Post.Prob Star
## 1      0.84     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

summary_all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  is_reactance,
  
  stats_to_report = c('CI'),
  model_rows_random = model_rows_random,
  model_rows_fixed = model_rows_fixed,
  model_rownames_random = model_rownames_random,
  model_rownames_fixed = model_rownames_fixed
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “is_reactance”

summary_all_models <- summary_all_models %>%
  print_df(rows_to_pack = rows_to_pack) %>%
  add_header_above(
    c(
      " ", "Hurdle Component" = 2, "Non-Zero Component" = 2,
      " " = 6
    )
  ) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 4,  
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      #"Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2
    )
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack,
  file.path("Output", paste0("AllModels", suffix, ".xlsx")), 
  merge_option = 'both', 
  simplify_2nd_row = FALSE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
summary_all_models
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Dichotome
Hurdle Component
Non-Zero Component
exp(Est.)_hu pa_sub 95% CI_hu pa_sub exp(Est.)_nonzero pa_sub 95% CI_nonzero pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log Est. mood_gauss 95% CI mood_gauss OR is_reactance 95% CI is_reactance
Intercept 2.49*** [ 1.53, 4.04] 49.28*** [41.95, 57.78] 123.18*** [109.03, 139.10] 3.80*** [ 3.58, 4.02] 0.57 [0.26, 1.24]
Within-Person Effects
Daily individual’s experienced persuasion 1.47*** [ 1.23, 1.79] 1.04 [ 0.99, 1.10] 1.02 [ 0.99, 1.06] 0.01 [-0.04, 0.06] 0.69* [0.48, 0.92]
Daily partner’s experienced persuasion 1.32** [ 1.10, 1.63] 1.02 [ 0.98, 1.07] 1.03 [ 0.99, 1.07] 0.01 [-0.04, 0.06] 1.12 [0.79, 1.66]
Daily individual’s experienced pressure 0.80 [ 0.47, 1.47] 0.89 [ 0.78, 1.00] 0.95 [ 0.86, 1.04] -0.05 [-0.18, 0.08] 2.31* [1.03, 9.08]
Daily partner’s experienced pressure 1.49 [ 0.81, 4.82] 0.95 [ 0.86, 1.05] 0.97 [ 0.90, 1.04] -0.06 [-0.19, 0.06] 1.79 [0.44, 13.34]
Daily individual’s experienced pushing 1.16 [ 0.87, 1.57] 0.99 [ 0.93, 1.05] 1.03 [ 0.97, 1.08] 0.04 [-0.04, 0.11] 1.43* [1.06, 2.07]
Daily partner’s experienced pushing 1.58** [ 1.17, 2.45] 0.98 [ 0.92, 1.04] 1.01 [ 0.96, 1.06] 0.11* [ 0.02, 0.18] 0.88 [0.54, 1.50]
Day 0.97 [ 0.63, 1.49] 1.01 [ 0.89, 1.16] 0.97 [ 0.88, 1.07] 0.24** [ 0.09, 0.39] 1.29 [0.45, 3.69]
Own actionplan NA NA NA NA NA NA NA NA NA NA
Partner actionplan NA NA NA NA NA NA NA NA NA NA
Daily barriers 1.21*** [ 1.10, 1.33] 1.04** [ 1.01, 1.07] 1.03** [ 1.01, 1.05] 0.13*** [ 0.10, 0.16] 1.03 [0.83, 1.29]
Daily weartime NA NA NA NA 1.00*** [ 1.00, 1.00] NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.87 [ 0.30, 2.51] 1.06 [ 0.77, 1.47] 1.02 [ 0.74, 1.39] 0.14 [-0.44, 0.72] 3.61 [0.51, 30.37]
Mean partner’s experienced persuasion 1.26 [ 0.43, 3.60] 0.98 [ 0.71, 1.34] 0.87 [ 0.63, 1.19] 0.43 [-0.15, 1.00] 6.05 [0.74, 67.12]
Mean individual’s experienced pressure 0.30 [ 0.06, 1.19] 0.66 [ 0.30, 1.47] 1.05 [ 0.74, 1.49] -0.18 [-0.78, 0.41] 42.14* [1.43, 2281.97]
Mean partner’s experienced pressure 0.25 [ 0.05, 1.01] 0.52 [ 0.23, 1.15] 1.04 [ 0.74, 1.47] -0.31 [-0.90, 0.28] 0.65 [0.01, 35.77]
Mean individual’s experienced pushing 1.13 [ 0.23, 5.53] 1.55 [ 0.86, 2.77] 1.12 [ 0.71, 1.78] 0.28 [-0.54, 1.12] 0.48 [0.02, 12.72]
Mean partner’s experienced pushing 1.49 [ 0.30, 7.30] 1.57 [ 0.87, 2.81] 1.29 [ 0.83, 2.05] -0.07 [-0.89, 0.77] 0.03 [0.00, 1.03]
Mean barriers 1.35* [ 1.07, 1.72] 1.08* [ 1.01, 1.16] 1.05 [ 0.99, 1.10] 0.36*** [ 0.28, 0.44] 1.10 [0.54, 2.40]
Mean weartime NA NA NA NA 1.00* [ 1.00, 1.00] NA NA NA NA
Random Effects
sd(Intercept) 1.20 [0.91, 1.63] 0.30 [0.22, 0.40] 0.32 [0.25, 0.43] 0.61 [0.48, 0.80] 1.95 [1.25, 2.94]
sd(Daily individual’s experienced persuasion) 0.16 [0.01, 0.44] 0.11 [0.06, 0.17] 0.05 [0.02, 0.09] 0.04 [0.00, 0.12] 0.47 [0.06, 0.93]
sd(Daily partner’s experienced persuasion) 0.22 [0.02, 0.51] 0.06 [0.01, 0.12] 0.05 [0.00, 0.10] 0.04 [0.00, 0.11] 0.40 [0.03, 0.99]
sd(Daily individual’s experienced pressure) 0.43 [0.02, 1.79] 0.07 [0.00, 0.26] 0.06 [0.00, 0.23] 0.07 [0.00, 0.27] 1.20 [0.12, 2.93]
sd(Daily partner’s experienced pressure) 0.77 [0.04, 2.34] 0.06 [0.00, 0.19] 0.04 [0.00, 0.13] 0.09 [0.00, 0.28] 1.64 [0.15, 4.12]
sd(Daily individual’s experienced pushing) 0.42 [0.04, 0.94] 0.05 [0.00, 0.14] 0.09 [0.02, 0.16] 0.07 [0.00, 0.18] 0.31 [0.02, 0.86]
sd(Daily partner’s experienced pushing) 0.42 [0.03, 1.15] 0.05 [0.00, 0.14] 0.05 [0.00, 0.12] 0.10 [0.01, 0.21] 0.42 [0.02, 1.22]
Additional Parameters
sigma NA NA 0.67 [0.65, 0.70] 0.54 [0.52, 0.56] 0.88 [0.85, 0.91] NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)

report::report_packages()
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • posterior (version 1.6.0; Bürkner P et al., 2024)
  • extrafont (version 0.19; Chang W, 2023)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • patchwork (version 1.3.0; Pedersen T, 2024)
  • R (version 4.4.2; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • ggridges (version 0.5.6; Wilke C, 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()